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I formulate a local density approximation for fermion sys- 
tems with pairing correlations based on a rapidly converging 
renormalization scheme for the pairing field. 

PACS numbers: 21.60. Jz, 21.30.Fe, 71. 15. Mb 



The proof that for any fermion system there ex- 
ist a unique energy density functional of the den- 
sity matter distribution p(r) alone, namely E gs — 
min J d 3 r£DFT(/?('")), is disarmingly simple How- 
ever, except for some trivial cases, the exact form of 
this functional is still a mystery and no constructive 
algorithms for its determination have been suggested 
so far. Significant progress has been achieved however 
within the Kohn-Sham Local Density Approximation 
(LDA) of the Density Functional Theory (DFT). For 
a normal fermion system (with no pairing correlations) 
Kohn and Sham have shown that the ground state en- 
ergy of any fermion system is a functional of its ki- 
netic energy and matter density distributions, namely 
E gs = min J d 3 r£LDA(7"(f ), p(r)). The current phi- 
losophy is that one should determine this functional 
from homogeneous infinite matter calculations and then 
use it to describe properties of either infinite inhomo- 
geneous or finite systems Q. By the same token, one 
would expect that the formulation of a LDA for fermion 
systems with pairing correlations should be straightfor- 
ward and that a corresponding universal LDA energy 
density functional £lda{t, p, v) of the kinetic energy 
r(r) = 2J2i |Vwi(r)| 2 , normal p(r) = 2£V k( r )| 2 and 
anomalous v{r) = v*{r)m{r) densities exists. (I shall 
be concerned here explicitly with the case of s-pairing 
in the so called weak coupling limit. Generalizations 
seem possible however.) The LDA extension described 
in Refs. is in terms of the anomalous density ma- 
trix v{ri,r<z) = (gs\ip^{r{)^pi{r2)\gs). Upon variation 
of the quasi-particle wave functions Vi{r),m{r) under 
standard restrictions one obtains the Kohn-Sham equa- 
tions, with a structure identical to the Hartree-Fock- 
Bogoliubov (HFB) or Bogoliubov-de Genes (BdG) equa- 
tions: 

[h{r) - p]u E {r) + A(r)v E (r) = Eu E (r), (1) 
A*{r)u E (r) - [h*{r) - p]v E (r) = Ev E (r), (2) 

where h(r) is the single-particle hamiltonian, A(r) = 
— gprfj and p is the chemical potential. Each quasi- 
particle state could be characterized by additional quan- 



tum numbers besides the quasi-particle energy E, which 
I shall not explicitly display however. In all the formulas 
presented here I shall likewise not display the spin degrees 
of freedom. One can show that the mere locality of the 
pairing field A(r) leads to a divergent diagonal part of 
the anomalous density v(r,r) HI). When |t*i — r%\ — » 
the anomalous density matrix has the singular behavior 
v{rx,r 2 ) = Y.E>o v *E( r 2) u E{ri) oc \ ri 1 _ r . 2 \- As a result, 
the local self-consistent pairing field A(r) cannot be de- 
fined. (When summing over the spectrum, the sum be- 
comes an integral if the spectrum is continuous and vice 
versa for an integral. I shall be casual in using either a 
summation or integration notation, hoping that the con- 
text makes this distinction obvious.) The existence of 
this particular divergence was the main obstacle in in- 
troducing an extension of the LDA approach to systems 
with pairing correlations. Fortunately, this divergence 
is one more example of the infinities which infest Quan- 
tum Field Theory (QFT) and for which the techniques to 
regularize them in a controlled fashion exist and can be 
extended and applied to inhomogeneous systems as well 
now. 

It is instructive to show how this divergence emerges 
and the simplest system to illustrate this is an infinite 
homogeneous one. Since the divergence is due to high 
momenta, thus small distances \r% — t*2 | , this type of di- 
vergence is universal and has the same character in both 
finite and infinite systems. Until recently methods to deal 
with this divergence were known only for infinite homo- 
geneous systems |6|-|l2| and only recently ideas were put 
forward on how to implement a renormalization scheme 
for the case of finite or inhomogeneous systems . As- 
suming for the sake of simplicity that the spectrum of 
the HF operator is simply e(k) = h 2 k 2 /2m, one can rep- 
resent the anomalous density matrix as follows || H 

d 3 k exp[ifc • (r% — t*2)]A 



v{r\,r 2 ) - 
d 3 k 



(27T) 



(2tt)3 2y /[ £ (k) p] 
3 exp[ik ■ (n - r%)] * — 



A 2 



(3) 



A 



A 



2[e(fc) - p - ij\ 



v 2y/[e{k) - p] 2 + A 2 
Amexp(z/ci?|ri — T2\) 
AirTi 2 \ri — r%\ 



(4) 



where p — h k F /2m. The last integral expression is well 
defined for all values of the coordinates r^a- Once one 
has recognized the existence of a divergence, the next 
step is to devise a way to regularize the theory. In a 
nutshell, what one has to do is to subtract the divergent 



part 



A-kTi 2 \ri —T2 | 



from the rest in the limit \r\ — r-2\ 



0. 
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Formally one can justify this apparently rather arbitrary 
procedure, either by following the steps outlined typically 
in renormalizing the gap equation in infinite systems - by 
relating the divergent part with the scattering amplitude 
[p|-|lO||- or by using well-known approaches in QFT - for 
example dimensional regularization [ jj"l|Jl^ ], or another 
QFT approach of introducing appropriate countcrtcrms 
with explicit cut-offs - or one can follow the philosophy 
of the pseudopotential approach ||,[||TJ] . In all cases one 
naturally arrives at the same final value for the gap. The 
renormalized gap equation can be written as 



1 

9 



d 3 k 
1 



1 



2y/[e(k)- p} 2 
ikpm 



2[e(fc) -p- i 7 ] 



where the coupling constant g is defined as 



gS(r 1 - r 2 ) = 



S 2 E, 



y* 



5v*(r l )8v(r 2 ) ' 



(5) 



(6) 



Previous approaches (6 12 use e(fc) only in the sec- 
ond term under the integral and in that case the last 
imaginary term does not appear. I have assumed 
here the simplest dependence of the LDA energy den- 
sity functional on the anomalous density v(r), namely 
£(r(r), p(r), \v(r)\ 2 ), merely for the sake of the simplic- 
ity of the presentation, but more general forms can be 
used as well. A note of caution: it would be incorrect to 
interpret some of the above formulas in the same manner 
as similar looking formulas appearing in various treat- 
ments of the pairing correlations with a zero-range in- 
teraction V(r\ — r 2 ) = gS(ri — r 2 ) (which can be related 
with the zero energy two-particle scattering amplitude 
g = 47rh 2 a/m). As it is well known for quite some time, 
even in the low density region, when kp\a\ -C 1, there are 
significant medium polarization corrections to the pair- 
ing gap The present LDA treatment is not limited 
by similar restrictions on the density. In the LDA energy 
density functional the polarization effects are already im- 
plicitly included in the definition of £(r(r), p(r), |i/(r)| 2 ) 
and the coupling constant g has no simple and direct re- 
lation to the vacuum two-particle scattering amplitude 
a. In this sense the LDA is similar in spirit to the Landau 
fermi liquid theory. 

Eq. (|^) can be used to extract from known proper- 
ties of homogeneous infinite matter (such as e(fe), A and 
density) the specific value of the coupling constant g to 
be used in constructing £(r(r), p(r), \v(r)\ 2 ). Assuming 
that a full microscopic calculation of homogeneous mat- 
ter at a given density p = kp /3ir 2 has been performed 
and that the value of the pairing gap at the fermi level 
is known, one can, using Eq. (^), calculate directly g(p) 
and thus obtain the simplest approximation to the LDA 
energy density functional £lzm(7~(t*), p( r ), W(r)\ 2 ) — 



£o(T~(r),p(r)) + g(p(r))\v(r)\ 2 , where £o(r(r), p{r)) is 
the Kohn-Sham energy density functional in the ab- 
sence of pairing correlations. In many treatments of 
the pairing correlations in infinite systems authors of- 
ten underline the dependence of the pairing gap on mo- 
mentum, that is A(fe). On one hand, typical calcu- 



lations 1 15 1 of the pairing field A(fe) in infinite sys- 
tems (with no medium polarization effects taken into 
account so far) show that for large momenta the pair- 
ing field decreases, as one would naturally expect. On 
the other hand, as soon as the momentum of a quasi- 
particle state is sufficiently different from the fermi 
momentum, when \k — kp\ » mA(kp)/h 2 kp <C kp, 
the effect of the pairing correlations on the single- 
particle properties is small, if not negligible. To a 
very good accuracy E(k) = (s(k) — p) 2 + A 2 (k) « 
\J (e(k) — p) 2 + A 2 (kp) « |e(fe) — p\ and thus the use 
of a fc-independent pairing field is a fair approximation. 
This is just another way of stating that the size of the 
Cooper pair h 2 kp/ mA || is much larger then the average 
interparticle separation in the weak coupling limit. Typ- 
ically this takes place when also the range of the pairing 
interaction is smaller than the size of the Cooper pair as 
well, and thus the pairing interaction could be described 
by a single coupling constant. 

Even though apparently the divergence has been suc- 
cessfully dealt with (in infinite homogeneous systems), 
a closer inspection of the entire approach reveals an in- 
consistency, which is somewhat hard to spot. The di- 
vergence is due to high momenta and for that reason 
one has subtracted the term A/2[e(k) — p — ij] in Eqs. 
(|],^]). Far away from the fermi surface however, the prob- 
lematic term A/2y / [e(fe) — p] 2 + A 2 behaves rather like 
A/2\e(k) — p\ instead. The main difference between these 
two subtraction procedures appears for hole-like states. 
As the fermi energy is finite, the integral over states below 
the fermi level is also finite. This feature, which breaks 
the approximate symmetry between the particle and hole 
states, is rather unsatisfactory and it has no theoretical 
underpinning. On one hand, in calculating the integral 
over the single-particle spectrum above the fermi level 
one expects a relatively fast convergence, when the en- 
ergy of the particle states is a "few gaps A away". On 
the other hand, the integral over the hole states con- 
verges only for energies of the order of the fermi energy 
ep = h 2 k F /2m. Clearly, in most cases of interest, the so 
called weak coupling limit, when A <C e_F, there is ab- 
solutely no physical reason to take into account single- 
particle states so far away from the fermi level in order to 
describe global or meanfield properties of nuclei in par- 
ticular. 

I show here how a relatively simple regularization 
scheme can easily deal with this problem in a very clear 
and easily implementable manner, suitable for any sys- 
tem, finite or infinite, homogeneous or inhomogeneous. 
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The regularized anomalous density is calculated from the 
following expression: 



* eg ( w 
H+E c 

fj,-E a 



dEg HFB {E)v E {r)u E {r) 



A(r)r E (r)Mr) Mr) 
2 |e-/i|+i7 2 



(7) 



r^(r- )M ), 



where gHFB(E) and gHF(s) are the HFB and HF density 
of states respectively, 



[h(r) - e]ip s (r) = 0, 
r(ri,r 2 ,^) = J deg H F(e) 



i)*(rx)i() s (r 2 ) 



(8) 
(9) 



2nh Iri — r 2 



s - /i| + 17 
r™»(r,/i) + 0(|ri-r 2 |). (10) 



7 is as usual a small infinitesimal quantity and r = T\% 
in the limit rq — r2| — > 0. As in Ref. I shall use 
a Thomas-Fermi approximation for the single-particle 
wave functions ip £ (r) and energies in order to evaluate 
the regulator. After introducing the local wave vectors 
l c {r) < k F {r) < k c (r) 



h(r) = - 


n 2 v 2 

2m + 


U(r), 


(11) 


lh 2 k 2 (r) 


+ U(r) 


= E c + fi, 


(12) 


2m 


\% 2 l 2 c {r) 


+ U{r) 


= -E c + fi, 


(13) 


2m 


h 2 k F (r) 
2m 


+ U(r) = 


= l l 


(14) 



and after some straightforward manipulations one can 
show that the renormalized anomalous density intro- 
duced above acquires the following form 



dEg HFB (E)v E (r)u E (r) 



A(r)mk c (r) 

2TT 2 h 2 
A(r)ml c (r) 

2n 2 h 2 



kjrjr) k c (r) + k F (r) 
' 2k c (r) n k c (r)-k F (r) 
k F (r) k F (r) + l c (r) 
2l c (r) U k F (r)-l c (r) 



(15) 



The only formal difference between this expression and 
the corresponding expression introduced in Ref. || is in 
the terms containing the second cut-off momentum l c (r) 
(last line). If either one of the wave vectors l c {r) or 
k c (r) becomes imaginary, then the corresponding terms 
in the renormalized anomalous density v reg (r) should be 
dropped. However, if the wave vector k F {r) becomes 
imaginary, the renormalized anomalous density is real 
and the above definition should be used, see Fig. |l] for a 
generic situation. 



It is convenient to introduce a notation for the cut-off 
anomalous density v c {r) :— J Q " dEgHFB{E)v* E (r)u E {r) 
and an effective position running coupling constant 



1 



9eff{r) 
ml c (r) 
2ir 2 h 2 



1 mk c (r 
9 



1 



k F {r) k c (r) + k F (r) 



2n 2 h 2 r 2fc c (r)~* k c {r) - k F [r) 
k F {r). k F (r)+l c (r) 



1 2l c (r) ln k F (r)-l c (r) 



(16) 



In the limit k F — > (in vacuum) the value of this effec- 
tive running coupling constant agrees with that derived 
in Ref. Jl^|. Using these notations one obtains for the 
renormalized pairing field 



A(r) = -gv r eg{r) = -g eff {r)v c (r). 



(17) 




FIG. 1. In region I all three wave vectors lc{r), k F (r) 
and k c (r) are real and both subtraction terms are present. In 
region II l c (r) is imaginary and the corresponding subtraction 
term in Eq. (^|) should be dropped. In region III all three 
wave vectors are imaginary and both subtraction terms in 
Eq. ( [l5| ) should be dropped. Even though in region II k F (r) 
becomes imaginary for larger r, the corresponding subtraction 
term containing k F (r) is real everywhere and it should be 
retained. 

Even though the cut-off momenta k c (r) and Idj") and 
the cut-off quasiparticle energy E c explicitly appear in 
the definition of both the effective coupling constant and 
of the cut-off anomalous density, the gap A(r) is indeed 
cut-off independent, once the cut-off energy E c has been 
taken sufficiently far from the fermi surface. This sit- 
uation is similar to the situation described in Ref. [gj, 
with the single difference that in the present case the 
convergence is achieved for significantly smaller values 
of E c . As one can judge from Fig. ^|, the present reg- 
ularization scheme is indeed very fast converging, while 
the regularization scheme presented in Ref. Q converges 
as expected at energies of the order of the fermi energy 
€ F = h 2 k E /2m. At the same time, the traditional ap- 
proach based on a (^-function with cut-off energy E c |ITJ] 
(for which g e ff(r) = g/[l — gmk c (r)/2ir 2 h 2 ]) converges 
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extremely slowly, and even at E c = 1000 MeV is still 
about 20% off the converged value. 

When computing the total energy of such a system one 
has to be careful and evaluate 



E„ 



2m 



T c (r) - A(r)u c (r) 



d 3 r£ (0,p(r)), 



(18) 



where the kinetic energy density is evaluated as r c (r) := 
2 J " dEgHFB{E)\VvE{r)\ 2 ■ (I have assumed here the 
simplest dependence of £o(r,p) on r.) Only this com- 
bined expression, containing the trace of the kinetic en- 
ergy with the trace of the pairing field and of the cut-off 
anomalous density, is converging as a function of the cut- 
off quasiparticle energy E c |o|. The reason is that r c (r) 
diverges in a similar manner as v c {f) as a function of E c . 
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FIG. 2. The gap A and the effective coupling constant g e g 
as a function of the cut-off energy E c for three regularization 
schemes. The full lines correspond to calculations using Eqs. 
(0~ 0)- Circles correspond to the regularization scheme pre- 
sented in Ref. [5] (when only terms with k c are present). The 
pentagrams correspond to the vacuum regularization scheme 
[16]. The calculation was performed for homogeneous neutron 
matter with p = 0.08 fm -3 and g = -250 MeV • fm 3 . 

The formalism described here paves the way to a LDA 
to pairing in the spirit of the Kohn-Sham theory . 
One has simply to add to the usual LDA energy density 
functional a pairing term g(p(r))\v c (r)\ 2 with a density 
dependent "bare coupling constant g(p(r))" , extracted 
from homogeneous infinite matter calculations. For the 
descriptions of many systems (e.g. nuclei, fermionic 
atomic condensates, 3 He and neutron matter) a term lin- 
ear in |y c (r)| 2 will most likely suffice. However, as we 
already know from the Landau-Ginzburg theory, terms 
proportional to |^ c (r)| 4 might become relevant and in 
such a case the energy density functional should be gen- 
eralized appropriately. Irrespective of the specific func- 
tional dependence of the energy density functional on 
the anomalous density v c (f), the emerging Kohn-Sham 
equations will be local and the ultraviolet divergence in 



the pairing field will have exactly the same character as 
the one studied here and consequently, can be dealt with 
using the same approach. 
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